This script is part 1 of our analysis of the stimulus-response characteristics of the SPARS, NRS_NP, and CNRS_P. This script generates exploratory plots of the relationship between stimulus intensity and scale ratings.

The three scales measure the following ranges of somatic sensation:

  • CNRS_P: 0 (no pain) to 100 (worst pain you can imagine)

  • NRS_NP: 0 (no sensation) to 100 (pain)

  • SPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)

Because the the stimulus range was calibrated to the sensitivity of each participant (compared to the fixed range of intensities used in Trial A), all analyses use the rank order of the nine stimulus intensities each participant was exposed to rather than the absolute intensities of the stimuli used.

For convenience, experimental blocks have been labelled in the order they took place within each scale. Each scale was tested across four successive blocks of 27 stimuli each. The order of in which the scales were assessed was randomized.


Import and inspect data

# Import
data <- read_rds('./data-cleaned/SPARS_B.rds')

# Inspect
glimpse(data)
## Observations: 2,256
## Variables: 9
## $ PID             <chr> "ID01", "ID01", "ID01", "ID01", "ID01", "ID01"...
## $ scale           <chr> "SPARS", "SPARS", "SPARS", "SPARS", "SPARS", "...
## $ block_number    <int> 2, 2, 2, 4, 4, 4, 6, 6, 6, 8, 8, 8, 11, 11, 11...
## $ trial_number    <int> 9, 15, 23, 7, 20, 25, 9, 18, 22, 3, 17, 23, 3,...
## $ intensity       <dbl> 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25, 2.25...
## $ intensity_char  <chr> "2.25", "2.25", "2.25", "2.25", "2.25", "2.25"...
## $ intensity_rank  <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ rating          <int> -31, -20, -48, -48, -21, -23, -48, -45, -47, -...
## $ rating_positive <dbl> 19, 30, 2, 2, 29, 27, 2, 5, 3, 0, 1, 3, 50, 50...

Clean and transform data

############################################################
#                                                          #
#                          Clean                           #
#                                                          #
############################################################
data %<>%
  # Select required columns
  select(PID, scale, block_number, intensity, intensity_char, 
         intensity_rank, rating, rating_positive) 

############################################################
#                                                          #
#                Calculate 'Tukey trimean'                 #
#                                                          #
############################################################
# Define tri.mean function
tri.mean <- function(x) {
  # Calculate quantiles
  q1 <- quantile(x, probs = 0.25, na.rm = TRUE)[[1]]
  q2 <- median(x, na.rm = TRUE)
  q3 <- quantile(x, probs = 0.75, na.rm = TRUE)[[1]]
  # Calculate trimean
  tm <- (q2 + ((q1 + q3) / 2)) / 2
  # Convert to integer
  tm <- as.integer(round(tm))
  return(tm)
}

############################################################
#                                                          #
#                    Generate core data                    #
#                                                          #
############################################################
# Add sequential block numbers, tagged by scale
data %<>%
    # Group by ID and scale
    group_by(PID, scale) %>% 
    # Arrange blocks
    arrange(block_number) %>%
    # Add sequential 
    mutate(block_sequential = dense_rank(block_number),
           block_sequential = paste0(scale, ' - ', block_sequential)) %>%
    ungroup()

# Convert NRS_NP and CNRS_P to SPARS-equivalent ranges (+50 to 50)
## Please see section: Equivalent unit rating for details
data %<>%
    # Fix column format to do the maths
    mutate(rating = as.numeric(rating)) %>%
    # CNRS_P conversion(divide rating by 2)
    # NRS_NP conversion (-100 and divide then by 2)
    mutate(rating_equivalent = case_when(
        scale == 'CNRS_P' ~ rating / 2,
        scale == 'NRS_NP' ~ (rating - 100) / 2,
        TRUE ~ rating
    ))

# Calculate the participant trimean for each scale (using 'rating_equivalent')
data_tm <- data %>% 
  group_by(PID, scale, intensity_rank) %>%
  summarise(tri_mean = tri.mean(rating_equivalent)) %>%
  ungroup()

Exploratory plots

Raw intensity ratings

Data are shown on their original scales: CNRS_P: 0 (no pain) to 50 (worst pain you can imagine), NRS_NP: -50 (no sensation) to 0 (pain), and SPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine).

# Generate plots
plot_raw <- data %>%
  group_by(PID) %>%
  nest() %>%
  mutate(plots = map2(.x = data,
                      .y = PID,
                      ~ ggplot(data = .x) +
                          aes(x = intensity_rank,
                              y = rating,
                              colour = scale,
                              fill = scale) +
                          geom_smooth(se = FALSE) +
                          geom_point() +
                          scale_y_continuous(limits = c(-50, 100)) +
                          scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                          scale_fill_brewer(name = 'Scale',
                                            type = 'qual', 
                                            palette = 'Dark2') +
                          scale_colour_brewer(name = 'Scale',
                                              type = 'qual', 
                                              palette = 'Dark2') +
                          labs(title = paste0(.y, ': Raw participant-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS'),
                               subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 100 (worst pain you can imagine)\nNRS_NP: 0 (no sensation) to 100 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                               x = 'Stimulus intensity (rank)',
                               y = 'Intensity rating') +
                          geom_hline(yintercept = 0, linetype = 2) +
                          facet_wrap(~ block_sequential)))

# Print plots
walk(plot_raw$plots, ~ print(.x))

Equivalent units ratings

Raw scores for the CNRS_P and NRS_NP scales (both rated on 0 to 100 scales) were converted to SPARS equivalent ranges (-50 to +50). The scaling of the CNRS_P and NRS_NP were as follows:

  • CNRS_P: raw 0 to 100 scores were converted to a 0 to 50 range by dividing 2, such that after the scaling, 0 = no pain and 50 = worst pain you can imagine. This is the positive range of the SPARS.

  • NRS_NP: raw 0 to 100 scores were converted to a -50 to 0 range by subtracting 100, and then dividing 2, such that after the scaling, -50 = no sensation and 0 = pain. This is the negative range of the SPARS.

# Generate plots
plot_equi <- data %>%
  group_by(PID) %>%
  nest() %>%
  mutate(plots = map2(.x = data,
                      .y = PID,
                      ~ ggplot(data = .x) +
                          aes(x = intensity_rank,
                              y = rating_equivalent,
                              colour = scale,
                              fill = scale) +
                          geom_smooth(se = FALSE) +
                          geom_point() +
                          scale_y_continuous(limits = c(-50, 50)) +
                          scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                          scale_fill_brewer(name = 'Scale',
                                            type = 'qual', 
                                            palette = 'Dark2') +
                          scale_colour_brewer(name = 'Scale',
                                              type = 'qual', 
                                              palette = 'Dark2') +
                          labs(title = paste0(.y, ': Scaled participant-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS'),
                               subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                               x = 'Stimulus intensity (rank)',
                               y = 'Scaled intensity rating (-50 to 50)') +
                          geom_hline(yintercept = 0, linetype = 2) +
                          facet_wrap(~ block_sequential)))

# Print plots
walk(plot_equi$plots, ~ print(.x))

Tukey trimean plots

# Generate plot
plot_tm <- data_tm %>%
    # Nest
    group_by(PID) %>%
    nest() %>%
    # Plot
    mutate(plot = map2(.x = data,
                       .y = PID,
                       ~ ggplot(data = ..1) +
                           aes(x = intensity_rank,
                               y = tri_mean) +
                           geom_smooth(method = 'loess',
                                       se = FALSE,
                                       colour = '#666666') +
                           geom_point(shape = 21,
                                      size = 3,
                                      fill = 'orange') +
                           scale_fill_brewer(name = 'Scale',
                                             type = 'qual', 
                                             palette = 'Dark2') +
                           scale_colour_brewer(name = 'Scale',
                                               type = 'qual', 
                                               palette = 'Dark2') +
                           scale_x_continuous(breaks = c(1, 3, 5, 7, 9)) +
                           scale_y_continuous(limits = c(-50, 50)) +
                           labs(title = paste0(.y, ': Scaled participant-level stimulus-response rating Tukey trimeans on the CNRS_P, NRS_NP and SPARS'),
                                subtitle = 'Plots are conditioned on the three scales\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
                                x = 'Stimulus intensity (rank)',
                                y = 'Scaled intensity rating (-50 to 50)') +
                           geom_hline(yintercept = 0, linetype = 2) +
                           facet_wrap(~ scale, ncol = 3)))

plot_tm$plot[[2]]

# Print plots
walk(plot_tm$plot, ~ print(.x))

Summary plot (equivalent units)

# Generate plot
data %>%
  ggplot(data = .) +
    aes(x = factor(intensity_rank),
        y = rating_equivalent,
        colour = scale,
        fill = scale) +
    geom_boxplot(alpha = 0.6) +
    scale_fill_brewer(name = 'Scale',
                      type = 'qual', 
                      palette = 'Dark2') +
    scale_colour_brewer(name = 'Scale',
                        type = 'qual', 
                        palette = 'Dark2') +
    labs(title = 'Scaled group-level stimulus-response ratings on the CNRS_P, NRS_NP and SPARS',
         subtitle = 'Plots are conditioned on experimental block\nCNRS_P: 0 (no pain) to 50 (worst pain you can imagine)\nNRS_NP: -50 (no sensation) to 0 (pain)\nSPARS: -50 (no sensation), 0 (pain threshold), +50 (worst pain you can imagine)',
         x = 'Stimulus intensity (rank)',
         y = 'Scaled intensity rating (-50 to 50)') +
    geom_hline(yintercept = 0, linetype = 2)


Session information

sessionInfo()
## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.3
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_GB.UTF-8/en_GB.UTF-8/en_GB.UTF-8/C/en_GB.UTF-8/en_GB.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2       forcats_0.2.0      stringr_1.2.0     
##  [4] dplyr_0.7.4        purrr_0.2.4        readr_1.1.1       
##  [7] tidyr_0.8.0        tibble_1.4.2       ggplot2_2.2.1.9000
## [10] tidyverse_1.2.1    magrittr_1.5      
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_0.12.15       RColorBrewer_1.1-2 cellranger_1.1.0  
##  [4] pillar_1.1.0       compiler_3.4.3     plyr_1.8.4        
##  [7] bindr_0.1          tools_3.4.3        digest_0.6.15     
## [10] lubridate_1.7.1    jsonlite_1.5       evaluate_0.10.1   
## [13] nlme_3.1-131       gtable_0.2.0       lattice_0.20-35   
## [16] pkgconfig_2.0.1    rlang_0.1.6        psych_1.7.8       
## [19] cli_1.0.0          rstudioapi_0.7     yaml_2.1.16       
## [22] parallel_3.4.3     haven_1.1.1        xml2_1.2.0        
## [25] httr_1.3.1         knitr_1.19         hms_0.4.1         
## [28] tidyselect_0.2.3   rprojroot_1.3-2    grid_3.4.3        
## [31] glue_1.2.0         R6_2.2.2           readxl_1.0.0      
## [34] foreign_0.8-69     rmarkdown_1.8      modelr_0.1.1      
## [37] reshape2_1.4.3     backports_1.1.2    scales_0.5.0.9000 
## [40] htmltools_0.3.6    rvest_0.3.2        assertthat_0.2.0  
## [43] mnormt_1.5-5       colorspace_1.3-2   labeling_0.3      
## [46] stringi_1.1.6      lazyeval_0.2.1     munsell_0.4.3     
## [49] broom_0.4.3        crayon_1.3.4